Загрузите данные в датафрейм.
gmp <- read.table("https://raw.githubusercontent.com/SergeyMirvoda/MD-DA-2017/master/data/gmp.dat")
gmp$pop <- gmp$gmp/gmp$pcgmp
estimate.scaling.exponent <- function(a, y0=6611, response=gmp$pcgmp,
predictor = gmp$pop, maximum.iterations=100, deriv.step = 1/100,
step.scale = 1e-12, stopping.deriv = 1/100) {
mse <- function(a) { mean((response - y0*predictor^a)^2) }
for (iteration in 1:maximum.iterations) {
deriv <- (mse(a+deriv.step) - mse(a))/deriv.step
a <- a - step.scale*deriv
if (abs(deriv) <= stopping.deriv) { break() }
}
fit <- list(a=a,iterations=iteration,
converged=(iteration < maximum.iterations))
return(fit)
}
a<-0.15
expr <- estimate.scaling.exponent(a)
estimate.scaling.exponent(a)
## $a
## [1] 0.1211533
##
## $iterations
## [1] 58
##
## $converged
## [1] TRUE
С помошью полученного коэффициента постройте кривую (функция curve) зависимости
plot(pcgmp~pop, data=gmp, xlab="Население", log="xy",ylab="Доход на душу населения ($/человеко-год)", main="Метрополии США, 2006")
curve(6611*x**expr[[1]],add=T,col='red')
Удалите точку из набора исходных данных случайным образом, как изменилось статистическая оценка коэффициента a?
deleting_row <- function(a,del.row = 0){
if(del.row==1){
row.for.del <- sample(nrow(gmp),1)
K <- which(!(1:nrow(gmp)) %in% row.for.del)
gmp <<-gmp[K,]
}
gmp$predicted <- 6611*(gmp$pop)^estimate.scaling.exponent(a)[[1]]
message(paste("a:",estimate.scaling.exponent(a)[[1]]," iterations:", estimate.scaling.exponent(a)[[2]],"\n"))
lines(gmp$pop, gmp$predicted, col = "red")
#return(estimate.scaling.exponent(a)[[1]])
}
Запустите оценку несколько раз с разных стартовых точек. Как изменилось значение a?
o1<-c()
i=1
for (u in seq(0.20,0.4,by = 0.01)){
o1[i]<-estimate.scaling.exponent(u)[[1]]
i<-i+1
plot(pcgmp~pop, data=gmp, xlab="Население", log="xy",ylab="Доход на душу населения ($/человеко-год)", main="Метрополии США, 2006")
deleting_row(u)
}
## a: 0.121153292739596 iterations: 70
## a: 0.121153292739592 iterations: 78
## a: 0.121153292739593 iterations: 100
## a: -0.0158227947406229 iterations: 100
## a: -0.282180595453594 iterations: 100
## a: -0.475270609311027 iterations: 100
## a: -0.719596634464836 iterations: 100
## a: -1.04635641522945 iterations: 100
## a: -1.48417713043446 iterations: 100
## a: -2.06949635226438 iterations: 100
## a: -2.85063412085794 iterations: 2
## a: -3.89182201309161 iterations: 2
## a: -5.27849638056654 iterations: 2
## a: -7.12437014241724 iterations: 2
## a: -9.58089019392308 iterations: 2
## a: -12.8498946401176 iterations: 2
## a: -17.2005650956782 iterations: 2
## a: -22.9921483818708 iterations: 2
## a: -30.7044336419058 iterations: 2
## a: -40.9786620271726 iterations: 2
## a: -54.6724802138931 iterations: 2
plot(o1,main="Оценка а с разных стартовых точек")
for(i in 1:100){
#print(paste("Количество точек", nrow(gmp)))
plot(pcgmp~pop, data=gmp, xlab="Население", log="xy",ylab="Доход на душу населения ($/человеко-год)", main="Метрополии США, 2006")
deleting_row(0.12,1)
}
## a: 0.121114444808736 iterations: 54
## a: 0.121109598973552 iterations: 54
## a: 0.121051791424099 iterations: 54
## a: 0.121063122493425 iterations: 54
## a: 0.121073882368287 iterations: 54
## a: 0.120962639733348 iterations: 54
## a: 0.120997617266159 iterations: 54
## a: 0.120982488200649 iterations: 53
## a: 0.120958073079915 iterations: 54
## a: 0.120969302525322 iterations: 54
## a: 0.120906250499022 iterations: 54
## a: 0.120951431111156 iterations: 54
## a: 0.120853195998464 iterations: 54
## a: 0.120802341916282 iterations: 53
## a: 0.120733852796318 iterations: 53
## a: 0.1204034114657 iterations: 53
## a: 0.120399875214202 iterations: 53
## a: 0.120556015906425 iterations: 53
## a: 0.120564855239912 iterations: 53
## a: 0.120491635720752 iterations: 53
## a: 0.120509338571741 iterations: 53
## a: 0.120504228779825 iterations: 53
## a: 0.120504165847206 iterations: 53
## a: 0.120549051575398 iterations: 53
## a: 0.120493073216661 iterations: 53
## a: 0.120614826927809 iterations: 53
## a: 0.120637954086077 iterations: 53
## a: 0.120611056710325 iterations: 53
## a: 0.120598440844703 iterations: 53
## a: 0.120537423988755 iterations: 53
## a: 0.12046843221171 iterations: 53
## a: 0.120509538771279 iterations: 53
## a: 0.120529594053472 iterations: 53
## a: 0.120508036439313 iterations: 52
## a: 0.120506086702969 iterations: 52
## a: 0.120499567611679 iterations: 52
## a: 0.120493177058957 iterations: 52
## a: 0.120479958735048 iterations: 52
## a: 0.120369508449169 iterations: 52
## a: 0.120386125745304 iterations: 52
## a: 0.120367624500289 iterations: 52
## a: 0.120349299709309 iterations: 52
## a: 0.120387027871802 iterations: 52
## a: 0.120372057746726 iterations: 52
## a: 0.12039504646698 iterations: 52
## a: 0.120407686223157 iterations: 52
## a: 0.120411112022311 iterations: 52
## a: 0.120341609610505 iterations: 51
## a: 0.120281050204919 iterations: 51
## a: 0.12032266622458 iterations: 51
## a: 0.120325902335481 iterations: 51
## a: 0.120379271312484 iterations: 52
## a: 0.120462893258827 iterations: 52
## a: 0.120468834620813 iterations: 52
## a: 0.120516915913398 iterations: 52
## a: 0.120460627673639 iterations: 52
## a: 0.120537524987198 iterations: 52
## a: 0.120516912270523 iterations: 52
## a: 0.12057491253572 iterations: 52
## a: 0.120612389670403 iterations: 52
## a: 0.120586085165935 iterations: 53
## a: 0.120506010814773 iterations: 52
## a: 0.120603503069724 iterations: 53
## a: 0.12057187863211 iterations: 53
## a: 0.120589879005083 iterations: 53
## a: 0.120620591644072 iterations: 53
## a: 0.120651081755946 iterations: 53
## a: 0.120624995371933 iterations: 53
## a: 0.120519862402356 iterations: 53
## a: 0.120558014926904 iterations: 53
## a: 0.120641067954741 iterations: 53
## a: 0.120903908632093 iterations: 54
## a: 0.120839183151217 iterations: 54
## a: 0.120884075745674 iterations: 54
## a: 0.120893377387248 iterations: 54
## a: 0.120881914533995 iterations: 53
## a: 0.120862005846508 iterations: 53
## a: 0.120848245022099 iterations: 53
## a: 0.120918111218219 iterations: 53
## a: 0.120845560403798 iterations: 53
## a: 0.120798376832012 iterations: 53
## a: 0.120638575858814 iterations: 53
## a: 0.120378650874897 iterations: 52
## a: 0.120442199792226 iterations: 53
## a: 0.120401481018942 iterations: 53
## a: 0.120388054749061 iterations: 53
## a: 0.120398773523536 iterations: 53
## a: 0.120349251907562 iterations: 52
## a: 0.120367062479791 iterations: 53
## a: 0.12039077103302 iterations: 53
## a: 0.12042560000916 iterations: 53
## a: 0.120317479135247 iterations: 52
## a: 0.120254098775309 iterations: 52
## a: 0.120209377634074 iterations: 51
## a: 0.120272847012106 iterations: 52
## a: 0.120234784940414 iterations: 51
## a: 0.120250921498287 iterations: 52
## a: 0.120285832190725 iterations: 52
## a: 0.12028358807766 iterations: 52
## a: 0.120349691629019 iterations: 53
а держится около 0.12 и количество итераций держится в одном интервале Количество итераций меняется при удалении значимых точек : чем ближе стартовое значение а к искомому, тем меньше количество итераций